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Numerical simulation using 
vorticity-vector potential formulation 

By H. Tokunaga 
1. Motivation and objectives 

An accurate and efficient computational method is needed for three-dimensional 
incompressible viscous flows in engineering applications. On solving the turbulent 
shear flows directly or using the subgrid scale model, it is indispensable to resolve the 
small scale fluid motions as well as the large scale motions. From this point of view, 
the pseudo-spectral method is used so far as the computational method. However, 
the finite difference or the finite element methods are widely applied for computing 
the flow with practical importance since these methods are easily applied to the 
flows with complex geometric configurations. However, there exist several problems 
on applying the finite difference method to direct and large eddy simulations. 

Accuracy is one of most important problems. This point has been already ad- 
dressed by the present author on the direct simulations on the instability of the 
plane Poiseuille flow and also on the transition to turbulence (Tokunaga, Ichinose 
& Satofuka, 1991a, b). In order to obtain high efficiency, the multi-grid Poisson 
solver is combined with the higher order accurate finite difference method (Toku- 
naga, Satofuka & Miyagawa, 1986). 

The formulation method is also one of the most important problems in apply- 
ing the finite difference method to the incompressible turbulent flows. The three- 
dimensional Navier-Stokes equations have been solved so far in the primitive vari- 
ables formulation. One of the major difficulties of this method is the rigorous sat- 
isfaction of the equation of continuity. In general, the staggered grid is used for the 
satisfaction of the solenoidal condition for the velocity field at the wall boundary. 
However, the velocity field satisfies the equation of continuity automatically in the 
vorticity-vector potential formulation. From this point of view, the vorticity-vector 
potential method was extended to the generalized coordinate system (Tokunaga, 
Yoyeda & Satofuka, 1991). In the present article, we adopt the vorticity-vector po- 
tential formulation, the generalized coordinate system, and the 4th-order accurate 
difference method as the computational method. At first, we present the computa- 
tional method and apply the present method to computations of flows in a square 
cavity at large Reynolds number in order to investigate its effectiveness. 

As is well known, the major drawback of the vorticity vector potential formulation 
is in the difficulty of specifying the boundary condition in the multiply connected 
domain. In applying the vorticity vector potential formulation to the flow with 
the complex geometric configuration, for example the flow along the multi -airfoil, 
we have to clear this hurdle. As the next step, therefore, we extend the present 
computational method to calculate the flow in a multiply connected domain. Lastly, 
the formulation of LES is dealt with in the framework of the present computational 
method. 
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2. Accomplishments 

2.1 Vorticity vector potential formulation 

Three-dimensional motions of fluid are governed by the Navier- Stokes equations 
and the equation of continuity 

du{ . . 1 A/ . x /i|X 

-Qjr +(U • vK = -VP+^ = z,y,*) (1) 

V # u = 0. (2) 

where u = (u x ,u y ,u z ) denotes the velocity, p the pressure, Re the Reynolds number, 
V the gradient operator, and A the Laplacian operator. Now we introduce the 
vorticity u) y the vector potential xj) y and the scalar potential <f> as 

u = v x 0 + V<t> ( 3 ) 

w = V x u (4) 

Then, the Laplace equation for the scalar potential, the vorticity transport equa- 
tions, and the Poisson equations for the vector potential are derived: 
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We introduce the generalized coordinate as 


* = *(£,»?), y = y(£,rj),z = *(C) 
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In order to specify the wall boundary condition, we define the the normal and 
the tangential component of ip as 
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tp n = n • %l>, 4>n = Ti • Vv, = r 2 • (!0) 

where n, n, and r 2 denote the normal and the tangential unit vector on the wall. 
Then, the wall boundary conditions for the vector potential are given as 

^=<( V , =^,=0 ( 1 ‘) 

On the other hand, the boundary vorticity is calculated from its definition 

u = v x u ( 12 ) 

2.2 Computational method 

The explicit method of lines is adopted as the computational method. In this 
method, spatial discretizations and the time integration are treated separately. For 
spatial discretization, the 2nd- and 4th-order accurate modified differential quadra- 
ture (MDQ) method is used. Partial derivatives of u;, for example, are approximated 
as 



L 
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The 2 nd-order accuracy is obtained by L = 1 and the 4th-order accuracy by L = 2 . 
Derivatives in the 77 - and (-direction are approximated in the same manner. After 
the above spatial discretization, the vorticity transport equations are reduced to a 
set of ordinary differential equations (ODEs) 

~ = F(u;) (15) 
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where /, , 7 , and I\ represent the grid number in the 77 -, and (-direction, respec- 
tively. As a time integration scheme of a set of ODEs, we apply the Runge-Kutta- 
Gill method. 

2.S Computation of flows in a square cavity 

In order to confirm effect ivity of the present computational method, we first carry 
out simulations of flows in a square cavity at Reynolds numbers Re = 10 3 and 10 4 . 
Figure 1 shows the velocity distribution on centerlines of the cavity. It is shown 
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FIGURE 1. Velocity distribution of the flow in a square cavity at Re = 10 3 and 
10 4 . Chain dot and chain dash line depict the result by the 4th- and 2nd-order 
method at Re = 10 3 respectively; solid and dash line represent the result by the 
4th- and 2nd-order method at Re = 10 4 respectively. Triangle and circle depict the 
results by Ghia et al. 

that the result of the 4th-order accurate method with the less grid points (65 x 65 
at Re = 10 3 ) is in good agreement with that of Ghia et al.(1982). 

In computation of the cavity flow at Re = 10 4 , we need the grid points 129 x 129 
in order to make the grid sufficiently dense in the vicinity of the wall. Figure 2 shows 
the convergence history of these computations with 2nd- and 4th-order accuracy. 
We needed 75000 time steps with At = 1/125 in order to obtain the converged 
solution when we used the 2nd-order accurate method. It is shown that the residual 



L<i -residual 
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oscillates drastically in time, and a great number of time steps are needed to attain 
the converged solution even for the 2nd-order method. However, the result with the 
4th-order accuracy shows that the L 2 residual preserves the constant level even at 
t = 900, and, therefore, it is concluded that the flow is unsteady. 



time step 

FIGURE 2. Convergence history of the cavity flow at Re = 10 4 . 

In order to investigate the temporal behavior of the flow, we show the vorticity 
contours at the early stage t = 50, 60, and 70, the middle stage t — 140, 150, 
and 160, and the later stage t = 240, 250, and 260. Since the present result has 
a fourth-order accuracy not only for the spatial discretization but also for time 
discretization, we can see the actual development of the cavity flow from the initial 
state in which the fluid is at rest. 

At an early stage, we find that there exists a number of vortices and that they are 
stretched in the course of time. The main vortex grows near the center of the cavity, 
and the other vortices Eire moved by rotation of the main flow. The stretching is 
caused by this main flow. However, the center of the main vortex is displaced by 
the mutual interaction of vortices. 

In the middle stage, we find that the main vortex is enclosed by the additional 
elongated vortex with the same sign. It is well known that a pair of vortices with the 
same sign rotate around each other along the weighting center of the pair. Thus, the 
unsteady motion of fluid is sustained. The center of the main vortex moves rapidly. 




t- 140 t = 150 t = 160 


FIGURE 4. Vorticity contours of the cavity flow at Re = 10 4 at middle time stage. 

The Biot-Savart velocity induced by the vortex is long-range, so the structure of 
small vortices in the corners of the cavity changes significantly in the course of time. 

The pair of vortices persist even at the later stage. However, at t — 260, the 
neighboring vortex is absorbed by the main vortex. Figure 6 shows the vorticity 
contours at t = 840 and 850. The vortex sheets are extremely stretched, and one of 
them is cut. The vortex generated by this mechanism interacts with other vortices, 
which explains the unsteadiness of thejlow in a cavity at high Reynolds number. 

2.4 Computation of flows in a multiply connected domain 

In general, the stream function value on the multiply connected domain is not 
known a priori. In order to resolve this problem, we apply the constraint proposed 
by Girault and Rivart (1979) 
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Figure 5. Vorticity contours of the cavity flow at Re — 10 4 at later time stage. 



t = 840 



t = 850 


FIGURE 6. Vorticity contours of the cavity flow at Re = 10 4 at the latest time 
stage. 

/f^ = ° (1?) 

r 

where T denotes the surface of the body placed in the flow. In the present compu- 
tation, the integral constraint is applied when the Poisson equation on the stream 
function is solved. We depict the stream line in Figure 7. The stream function value 
is specified to 1 on the upper wall and 0 on the lower wall. We ultimately obtained 
the stream function value -0.11375 on the square, which shows a good agreement 
with the result of Lippke and Wagner (1091). 
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FIGURE 7. Streamline lines of the flow in a multiply connected domain at Re = 20. 


2.5 Computation of transition of channel flow using LES 

For the subgrid-scale model, we choose the dynamical model proposed by Ger- 
mano et al (1991). Then, the basic equation of LES is obtained in the vorticity 
vector potential formulation as 
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where Vt denotes the dynamical subgrid scale turbulent viscosity and Sij the strain 
tensor. The generalized coordinate transformation is applied to this equation. 

At first we will carry out LES of a transition of plane channel flow in order to 
investigate the validity of the present method. The initial condition is created by 
using the result of the direct simulation of 2-D and 3-D linear instability (Tokunaga, 
Satofuka & Miyagawa, 1986) 


3. Summary and future plans 

An accurate and efficient computational method is presented for incompressible 
viscous flows. It is shown that the present method well predicts the flow in multiply 
connected domain. The 4th-order accurate method shows that the square cavity 
flow can be calculated accurately by less grid points in comparison with the 2nd- 
order method at Re = 10 3 . Further, it is shown that the cavity flow is unsteady 
at Re = 10 4 , and actual process of the flow development is cleared by using the 
4th-order accurate method. 

The next step in this work will be to test the large eddy simulation of transition 
in a plane channel dealt with by Germano et al. (1991), and this study is under 
way. The ultimate goal of this work is direct and large eddy simulation of the flows 
in a multiply connected domain, which is of practical importance in aerodynamics. 
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